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We develop a Green's function approach to quasiparticle excitations of open-shell systems within 
the GW approximation. It is shown that accurate calculations of the characteristic multiplet struc- 
ture require a precise knowledge of the self energy and, in particular, its poles. We achieve this by 
constructing the self energy from appropriately chosen mean-field theories on a fine frequency grid. 
We apply our method to a two-site Hubbard model, several molecules and the negatively charged 
nitrogen-vacancy defect in diamond, and obtain good agreement with experiment and other high- 
level theories. 



Introduction. — In nature, there exists a wide range of 
electronie systems with open shells, including most atoms 
and many molecules, but also defects in crystalline solids. 
These systems play important roles in almost all areas of 
physics, chemistry and biology: for example, the neg- 
atively charged nitrogen- vacancy (NV~) defects in dia- 
mond are used for biological imaging [l|, Q and are also 
promising candidates for qubits in quantum computers 
[3-5]. 

It is therefore important to develop theoretical meth- 
ods to study open-shell systems and their properties. 
While for closed-shell systems a well-established set of 
methods exists, ranging from wave function-based quan- 
tum chemistry approaches to density-functional theory 
(DFT) and Green's function based many-body perturba- 
tion theory, the accuracy of these methods when applied 
to open-shell systems is less certain: Even the applica- 
tion of wave function-based methods to small open-shell 
molecules is far from straightforward 0] and standard 
density functionals are known to break the orbital and 
spin degeneracy of the ground state 0, S] ■ 

A Green's function approach to electron excitations in 
open-shell system was first considered by Cederbaum and 
coworkers [3, [13] in the 1970's. These authors only ap- 
plied the formalism to toy problems with few orbitals and 
employed approximations to the self energy which are 
not feasible for large systems, such as open-shell defects. 
Other previous applications of Green's function theory 
to open-shell systems using the GW approximation had 
either carefully selected reference states to avoid compli- 
cations associated with the open-shell (ill or ignored the 
degenerate ground-state problem [12!, Il3|. 

In this Letter, we extend the GW approach to open- 
shell systems. Calculations on several prototypical sys- 
tems are performed: a two-site Hubbard cluster, four 
molecules (nitrogen dioxide, oxygen, nitrogen difluoride, 
chlorine dioxide) and the NV~ center in diamond. We 
find our approach is capable of describing these systems 
with quantitative accuracy. We have identified and im- 
plemented two important elements for accurate results in 
GW calculations of open-shell systems: i) a careful choice 



of the mean-field starting point providing accurate self- 
energy pole positions, and ii) a method for evaluating the 
self energy on a fine frequency grid. 

Theory. — In a photoemission experiment with photons 
of energy Wphoton (setting h — 1), the photocurrent J(efe) 
due to photoelectrons with momentum k and energy ek 
is given by [l3] 



J(,(-k) — ^ AfciAjfc^y (efc — Wphoton), 



(1) 



where A^^ = (fe] AdipoidV'j) and Ay(w) = 
{ipi\A{r,r',uj)\ipj) denote matrix elements of the 
dipole operator and the spectral function, respectively, 
with -tpi being an appropriate single-particle orbital. 
Neglecting off-diagonal matrix elements for an ap- 
propriately chosen physical set of orbitals, we obtain 
Ajj{Lu) = l/7r]ImG'jj(a;)] by computing the interacting 
Green's function (here we give the electron removal 
part) 



G«(^) = E 
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with Ex = E^g^^ - e[^ ^\ Here, \N, 0) and E^^^^ denote 
the A^-particle ground state and its energy, respectively, 
while \N — 1, A) denotes an {N — l)-particle state (with 
A being an appropriate set of quantum numbers) with 
energy e'^ ^\ Also, Cj is the destruction operator for 
an electron in orbital j and ?/ = 0+. 
E\ solves the quasiparticle equation 



Ex = e. 



^,,iEx)-V^, 



(3) 



where ej and Vjj^ 



denote the orbital energy and a diago- 
nal matrix element of the exchange-correlation potential 
from a mean-field calculation, respectively, while (cj) 
is a diagonal matrix element of the self-energy operator. 

The quasiparticle equation jEq. ([3|)] follows from 
Dyson's equation [H] 



(4) 
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which relates the interacting Green's function to the 
mean-field Green's function G'o,y(w) via the self energy. 
The standard derivation of Dyson's equation [31 assumes 
the existence a nondegenerate interacting ground state 
which evolves into a nondegenerate single Slater deter- 
minant state as the interactions are adiabatically turned 
off. The hallmark of open-shell systems, however, is the 
existence of multiple degenerate ground states which do 
not generally evolve into noninteracting single Slater de- 
terminant states [3| ■ If — for a particular ground state 
— the resulting noninteracting state is a sum of Slater 
determinants, one has to employ the methods of quantum 
field theory with initial correlations and replace Dyson's 
equation with a more complicated expression [la . [l7| . 
In our calculations, we avoid this difficulty by carefully 
choosing a ground state which evolves into a single Slater 
determinant such that Dyson's equation is valid. In par- 
ticular, we work with the ground state with the highest 
magnetic quantum number because there exists a corre- 
sponding single Slater determinant with the same prop- 
erties (i.e., it is also an eigenstate of the total spin and/or 
orbital angular momentum operator with the same eigen- 
value) jlOl] . An approximation to this particular ground 
state is provided by standard spin-polarized mean-field 
calculations. We note that it is not always possible to 
find a single determinant ground state. However, such a 
state must exist whenever Hund's rules apply. 

In closed-shell systems, Eq. ([3]) typically has a single 
solution leading to a pronounced quasiparticle peak in 
Ajj{ui) which corresponds to the removal of an electron 



from orbital j [18|. In open-shell systems, the orbital and 
spin angular momenta of the electrons in the unfilled 
shells can couple in various ways resulting in multiple 
low-energy eigenstates of the N and the {N — l)-particle 
system. The coupling of angular momenta generally pro- 
duces eigenstates which are sums of multiple Slater deter- 
minants [l9| . As a consequence, multiple eigenstates of 
the (A'^— l)-particle system can make significant contribu- 
tions to Gjj (w) if their matrix element in the numerator 
of Eq. ([2]) is large. Gjj (lu) then has multiple poles and we 
expect to find multiple solutions of Eq. ([3]) . This impor- 
tant connection between the poles of the self energy and 
the multiplet structure of open-shell systems was first 
established by Cederbaum and coworkers [^, ■ 

If Gjj (w) has multiple poles, Eq. (jU shows that the self 
energy Sjj (w) must also have poles occurring between the 
poles of Gjj{uj). The occurrence of poles in Tijj{uj) near 
Ex is a particular feature of open-shell systems and a 
direct consequence of the electronic multiplet structure. 

In actual calculations for open-shell systems, a precise 
knowledge of the frequency dependence of the self energy 
is necessary to locate its poles and obtain accurate mul- 
tiplet splittings. In contrast, for closed-shell systems it 
is usually sufficient to employ a simple linear expression 
for the frequency dependence of the self energy in the 
vicinity of the quasiparticle energy [l8| . 



In this work, we employ the GW approximation to 
the self energy following the first-principles method of 
Hybertsen and Louie [l8|. To obtain Yijj{u)) at many 
frequencies, we make use of a specific form of the eval- 
uation of the frequency dependence of the dielectric re- 
sponse and self energy as proposed in Refs. and 11 1. 



In this approach, Ejj(w) is separated into a frequency- 
independent bare exchange part Ej^"* and a frequency- 
dependent correlation part Ej^-^(w) given by 



jnl I 



nl 



(5) 



where /i denotes the chemical potential and f2/ is a neu- 
tral excitation energy of the A^-particle system obtained 
by solving Casida's equation in the random-phase ap- 
proximation [20| . Also, Vjni denotes a Coulomb matrix 
element between t he p roduct '4'j''Pn and the fluctuation 
charge density pi [20| (see Supplementary Material for 
details on the approach). 

Equation ([5]) shows that the poles of Ejj(aj) are deter- 
mined by the mean-field electron removal (or addition) 
energies e„, which are the poles of Gq, and by the neu- 
tral excitation energies fi/, which are the poles of the 
screened interaction Wq in the random-phase approxima- 
tion. Both Era and flj depend on the mean-field theory 
used to compute Go and Wq, implying an analogous de- 
pendence on the choice of the mean-field starting point 
for the poles of 'Ejj{uj). 

In principle, the self energy should be computed from 
the interacting Green's function G, whose poles are at 
Ex, and the exact screened interaction W For 
closed-shell systems, it is possible to carry out self- 
consistent GWo calculations where the self energy is re- 
computed using the iterated Green's functions such that 
E becomes independent of the mean-field starting point 
21| . For open-shell systems, self-consistent calculations 
are more difficult because of the more complicated struc- 
ture of G and additional problems to be discussed below. 
To obtain accurate self-energy pole positions we instead 
carefully choose mean-field theories that yield e„ and fi/ 
which are good approximations to Ex and the poles of the 
exact W, respectively. In general, one finds that the poles 
of Wq obtained from standard density-functional calcula- 
tions are good approximations to neutral excitation ener- 
gies. In contrast, the poles of Gq obtained from density- 
functional theory often differ from the exact removal or 
addition energies (i.e. the quasiparticle energies) by sev- 
eral electron volts. Such an error in the poles of Go leads 
to a similar-sized error in the self energy pole locations 
and to a large error in the multiplet splittings. To obtain 
the best Gq, we construct it from mean-field calculations 
using the static COHSEX approximation 22, 23l |. 

In addition, if the result of a calculation depends on 
a particular self-energy pole we carry out partially self- 
consistent calculations where we only update the partic- 
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ular En in Eq. ([5]) which determines the position of the 
self-energy pole under consideration. 

Molecules. — First, we study the electronic multiplet 
structure of four small molecules for which accurate ex- 
perimental data is available. 

Nitrogen dixoide (NO2) has a doublet ground state. 

251 at the ex- 



TABLE I: Comparison of our results for NO2 with experiment 
[2^ . All energies are given in eV. 
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We first carry out DFT calculations 
perimental geometry [2^ using the spin-polarized LDA 
exchange-correlation functional, norm-conserving pseu- 
dopotentials, a plane-wave basis (50 Ry cutoff) and a 
cubic supercell with linear dimension of 10.6 A. 

For the construction of Wq we use wave functions and 
energies from the DFT calculation. We use 300 empty 
states and a 15 Ry momentum space cutoff for the dielec- 
tric response. For Go we use wave functions and energies 
from a static COHSEX calculation. Table U shows that 
the COHSEX single-particle energies are much closer to 
the experimental ionization potentials than the DFT en- 
ergies, but the multiplet structure is still missing in this 
calculation. For the calculation of the self-energy matrix 
element we use 300 empty states and a modified static 



remainder correction |27l . |28| which extends the sum over 
n in Eq. ([S]) to all empty states and greatly improves con- 
vergence. This choice of parameters results in multiplet 
splittings converged to within ~ 0.1 eV. 

Figure [TJa) shows the self energy and spectral func- 
tion for the removal of a down-spin electron from the 4&2 
orbital [see insert in Fig. [UJa)]. We do not expect any 
multiplet structure for this process because the up-spin 
hole can only couple to the up-spin electron in the 6a 1 
orbital to give a triplet state. Indeed, the spectral func- 
tion exhibits a single peak corresponding to the triplet 
(^Sa) state. 

Figure [Ijb) shows results for the removal of an up- 
spin electron from the 462 orbital. The down-spin hole 
can now couple to the up-spin electron in the 6ai or- 
bital to yield either a singlet (^i?2) or a triplet ('^-62) 
state. Indeed, we find two solutions of Eq. ^ resulting 
in two poles of the Green's function and two peaks in the 
spectral function with a singlet-triplet splitting of 1.8 eV 
which compares favorably with the experimental splitting 
of 1.5 eV (Table HI). In contrast, the singlet-triplet split- 
ting from Glda^lda is 2.9 eV highlighting the impor- 
tance of an accurate mean-field starting point. To make 
sure that the two solutions are indeed multiplet states we 
traced back the low lying self-energy pole to open-shell 
features in Go and Wo: namely, to the pole in Go due to 
the unpaired up-spin 6ai state and the pole in Wq due to 
the 4b2| — > 6bi4. transition between the two open shells. 
Hund's rule suggests that the lower energy solution is the 
triplet state. 

Inspection of Table U shows that we obtain two values 
for the energy of the triplet state '^62, one from the re- 
moval of an up-spin electron from the 4^2 orbital, one 
from the removal of a down-spin electron from the same 
orbital. These values differ by 0.8 eV and bracket the ex- 
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FIG. 1: Self energy {uj) and spectral function Ajj (fx?) for 
(a) the removal of a down-spin electron from the j — 462 
orbital in NO2 and (b) the removal of an up-spin electron 
from the j = 462 orbital. A Lorentzian broadening of 20 meV 
is used for each curve. 



perimental result. There are two factors which contribute 
to this discrepancy: i) remaining errors in the positions of 
the self-energy poles which contaminate only solutions of 
the up-spin quasiparticle equation and ii) missing vertex 
corrections which contaminate solutions of the up- and 
down-spin quasiparticle equations in different amounts 
[2^. We expect that the inclusion of vertex corrections 
will reduce the difference. Nevertheless, as shown above, 
accurate multiplet splittings can be extracted from our 
calculations if the energy differences are calculated from 
solutions of the quasiparticle equation for a particular 
spin direction. 

The ratio of the areas under the singlet and the triplet 
peaks in Fig. [IJb) should be the experimentally observed 
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TABLE II: Comparison of our results for O2, NF2 and CIO2 
with experiment [31I - I33I ]. All energies are given in eV. 
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ratio of photoemission intensities, the so-called multiplet 
ratio |30|. We find in our calculations that the multi- 
plet ratios are much more sensitive to the positions of 
the self-energy poles than the multiplet splittings. We 
do not expect that these ratios can be computed reliably 
with our current GW approach because of the remaining 
uncertainties in the self-energy pole locations. Hovifever, 
Schirmer and coworkers found a relatively simple analyt- 
ical procedure for calculating these ratios based on the 
addition of angular momenta [s^l- We expect that the 
combination of their approach for the multiplet ratios 
and the GW approach for the multiplet splittings offers 
a reliable and complete description of the multiplet struc- 
ture of open-shell systems. 

Table |lT] shows our results for the oxygen (O2), ni- 
trogen difluoride (NF2) and the chlorine dioxide (CIO2) 
molecules. The GW multiplet splittings are 2.4 eV for 
O2, 2.3 eV for NF2 and 2.5 eV for CIO2. They compare 
favorably with experimental splittings: 2.3 eV for O2, 
1.8 eV for NF2 and 2.4 eV for CIO2 [31-33]. However, 
splittings obtained from Glda^lda can deviate from 
experimental findings by several electron volts. 

Hubbard cluster. — To further establish the accuracy of 
our method, we apply it to an analytically solvable model 
system: a two-site Hubbard cluster containing three elec- 
trons with a Hilbert space spanned by four spin-orbitals. 
We denote the hopping parameter t and the on-site in- 
teraction U. This system has a doublet ground state. 
We compute the mean-field wavefunctions and energies 
using the spin-polarized Hartree-Fock method and then 
evaluate the self energy corresponding to the removal of 
an up-spin electron from the doubly occupied bonding 
orbital. We find two solutions of the quasiparticle equa- 
tion due to the occurence of a pole in the self energy: 
their separation is I.IU for U/t < 1 where we expect the 
GW approximation to give accurate results. 

In this model system, analytical evaluation of the 
Hamiltonian for the three and two particle systems al- 
lows for the calculation of the exact many-body Green's 
function which agrees well with the GW result: it has 
two poles corresponding to a singlet and a triplet state 
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FIG. 2: Self energy (cj) and spectral function A.jj (iu^ for 
the removal of an up-spin electron from the j = orbital. A 
Lorentzian broadening of 5 meV is used for each curve. 



TABLE III: Comparison of the calculated multiplet splittings 
for the NV~ defect in diamond with results from exact diag- 
onalization calculations on the extended Hubbard modelj35l]. 
All energies are given in eV. 
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of the two particle system separated by U. 

NV~ center. — Next, we apply our approach to the 
NV~ center in diamond which has a triplet ground state. 
This defect complex currently attracts much attention 
because of its extraordinary properties, such as long co- 
herence times and potential application to quantum com- 
puting 

Again, we carry out DFT calculations as described in 
the previous sections. We employ a 64-atom supercell 
and relax all ionic positions. As a test, we first carry out 
GW calculations with a generalized plasmon pole model 
[3^ ] and find good agreement with the similarly calcu- 
lated results of Ma and Rohlfing [l^j, who use a 256- 
atom supercell, for the energy differences of the defect 
levels in the gap. This indicates that our supercell size is 
sufficient to extract accurate multiplet splittings for the 
defect levels using the current method. 

For the NV~ center we find that the DFT orbital ener- 
gies are much closer to the static COHSEX results than 
in NO2, and we can use the DFT energies and wave func- 
tions both for the construction of Wq and Gq. We use 
300 empty orbitals and a 12 Ry momentum space cut- 
off for Wq and 300 empty orbitals in conjunction with a 
modified remainder correction for the self-energy matrix 
elements. This choice of parameters results in multiplet 
splittings converged to within ~ 0.1 eV. 

Figure [5] shows results for the removal of an up-spin 
electron from the i/ level. As in the NO2 calculation, the 
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self energy exhibits a low-lying pole leading to two solu- 
tions of the quasiparticle equation. To understand which 
many-body states these solutions correspond to, we com- 
pare our results to exact diagonalization calculations of 
the extended Hubbard model of Choi, Jain and Louie (35|. 
These authors fit the parameters of an extended Hubbard 
model for the defect levels to ab initio static COHSEX 
results and show that this model describes accurately 
neutral excitations. The model predicts four many-body 
states 5' A for the v^e'^ configuration. However, only two 
of the four states, namely and ^^2, are observed in 
our calculations because by symmetry only these states 
have a non- vanishing matrix element xlc^-^lv-fV^ex-i-eyi-) 
with the ground state. Table IHII shows that we obtain 
good agreement with the extended Hubbard model re- 
sults for the multiplet splittings. Note that the ^A2-^^2 
splitting in Fig. [2] corresponds to the last row in Table Hill 
The first row in Table Hill shows the splitting between the 
'^E state obtained by removing an up-spin electron from 
the e defect orbital and the state. Again, we find 
good agreement with the extended Hubbard model re- 
sults. 
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